Generating scanning spot locations for laser eye surgery

ABSTRACT

Scanning spot locations are generated for ablating tissue using a scanning laser beam over a treatment region by fitting a target function representing a desired lens profile of ablation with a basis function representing a treatment profile produced by overlapping scanning spots in a particular treatment pattern. Symmetry effects are utilized to simply the process for determining the scanning spot locations. In some embodiments, the basis function is a two-dimensional function representing a two-dimensional section of a three-dimensional treatment profile which has symmetry with respect to the two-dimensional section extending along the treatment pattern. For example, the treatment pattern is generally straight for myopic and hyperopic cylinders, and is generally circular for myopia and hyperopia. The target function and the basis function may be discrete for implementation in a software algorithm, and be fitted using a least square fit. The fit produces ablation depths for discrete scanning spots which are used to calculate the number of pulses at each reference position along the two-dimensional section. The pulses are distributed along the treatment pattern to produce the desired overlapping effect.

[0001] This application is based on and claims the benefit of U.S.Provisional Patent Application No. 60/189,633, filed Mar. 14, 2000, theentire disclosure of which is incorporated herein by reference.

BACKGROUND OF THE INVENTION

[0002] The present invention relates generally to tissue removaltechniques and, more particularly, to generating locations for scanningwith a laser to achieve a desired ablation profile for correction oferrors in vision during laser eye surgery.

[0003] A scanning system has the ability to trace out an arbitrarypattern with a small low energy spot. In most cases, a small spotequates to finer scanning details but at the expense of requiring morepulses to remove a given volume. The notion that a small spot will givea better fit than a large spot is generally true for arbitrary spot andablation shapes, but the spot shapes can also affect the fit. Forinstance, trying to fit round disks into a square shape will result in aresidue. The emphasis on getting a good fit should be on choosing a goodbalance of spot geometry and size.

[0004] Small spot scanners have their problems. A small spot will havesmaller coverage per shot, thereby requiring more pulses to remove agiven ablation volume (inversely proportion to the area of the spotsize). A larger spot will require a substantially smaller number ofpulses but at the expense of resolution. Understanding how treatmentvaries with spot overlap will ease our ability to create ablationpatterns.

[0005] Present ablation algorithms follow a removal or subtractiveprocess. An ablation spot is placed on a location along thetwo-dimensional corneal surface. This spot ablates a volume of tissuefrom the surface to produce a crater on this surface. Another spot isapplied at another location along the two-dimensional corneal surfaceand another crater is produced. The challenge lies in placing thecraters in the correct locations such that adding up the overlappingcraters will produce the desired surface without producing undesirableresidues and requiring excessive processing time.

SUMMARY OF THE INVENTION

[0006] The present invention relates to tissue ablation utilized in, forinstance, corneal sculpting. In general, embodiments of the inventiondetermine a treatment table containing scanning spot locations andcharacteristics (e.g., size, shape, and depth) for directing overlappingscanning spots of a laser beam to achieve a target ablation profile.Some embodiments of the invention use symmetry effects (e.g.,axisymmetry or bilateral symmetry) to simplify the overlap analysis todetermine scanning spot locations for directing the laser beam toachieve a desired ablation profile. In specific embodiments, the spotcenter locations define a series of linear overlapping ablation paths,with the shape of the ablation paths preferably being selected inresponse to a desired change in an optical characteristic of the eye.The scanning spots can be uniform or variable in size, shape, and/ordepth. Advantageously, these paths can produce ablation treatmentprofile elements which can be combined, adjusted, and positioned overthe two-dimensional corneal surface to, in combination, produce thedesire three-dimensional sculpting. The ablation profiles of theoverlapping paths can be represented by basis functions, which mayinclude an array of ablation profile data or can be expressedanalytically for certain profiles.

[0007] An aspect of the present invention is directed to a method ofgenerating a treatment table for ablating tissue using a scanning laserbeam for generating scanning spots over a treatment region larger inarea than the scanning spots. The method comprises providing a targetfunction or lens function representing a desired lens profile forablating the tissue by scanning spots of the laser beam on the tissue. Abasis function represents a treatment profile produced by theoverlapping scanning spots in a treatment pattern. The target functionis fitted with the basis function to obtain the treatment tableincluding scanning spot locations and characteristics for theoverlapping scanning spots of the laser beam. In a specific embodiment,the scanning spot locations are randomized to produce a random scanningorder. The scanning spot characteristics of a scanning spot may includescanning spot size, shape, and depth of the scanning spot at a specificscanning spot location.

[0008] Any desirable treatment pattern for scanning with overlappingscanning spots of the laser beam may be specified. In some embodiments,the overlapping scanning spots define a treatment pattern in the formof, e.g., a linear path, and the basis function is a two-dimensionalfunction representing a two-dimensional section of a three-dimensionaltreatment profile which has symmetry with respect to the two-dimensionalsection extending along the treatment pattern. For example, thetreatment pattern is generally straight for myopic and hyperopiccylinders, and is generally circular for myopia and hyperopia. In aspecific embodiment, the lens function represents an ablation depth as afunction of a distance from an optical axis of a cornea. The basisfunctions of a series of offset patterns are combined, the profiles oflaterally adjacent paths often overlapping to provide a smooth treatedsurface.

[0009] In some embodiments, the spot size and shape are generally fixed.In other embodiments, the spot size and shape are variable. The abilityto achieve the desired ablation profile using fixed spot size and/orshape may allow the use of a more simplistic laser source and simplifythe ablation process.

[0010] In specific embodiments, fitting the lens function with the basisfunction includes fitting at N discrete evaluation points. The basisfunction includes M discrete basis functions representing M overlappingscanning spots. The M discrete basis functions may represent Moverlapping scanning spots across a treatment zone length representingthe length across a generally two-dimensional section which is orientednormal across a generally straight treatment pattern or which isoriented radially across a generally circular treatment pattern.

[0011] For a treatment profile having a generally uniformtwo-dimensional section oriented normal across a generally straighttreatment pattern, the discrete basis functions represent thetwo-dimensional section as

X _(i)(x _(j))=y _(i)(x _(j))={square root}{square root over ((s/2)²−(x_(j) −x _(0i))²)}; or,

[0012] for a treatment profile having a generally uniformtwo-dimensional section oriented radially across a generally circulartreatment pattern, the discrete basis functions represent thetwo-dimensional section as${X_{i}\left( x_{j} \right)} = {{\theta_{i}\left( x_{j} \right)} = {\cos^{- 1}\left( \frac{x_{j}^{2} + x_{0i}^{2} - \left( {s/2} \right)^{2}}{2 \cdot x_{0i} \cdot x_{j}} \right)}}$

[0013] where

[0014] s is the diameter of the scanning spot;

[0015] j=1, . . . , N;

[0016] x_(j) is a reference x-coordinate for the two-dimensional sectionmeasured from an optical axis of the cornea of a j^(th) evaluation pointfor the center of the scanning spot;

[0017] x_(0i) is an x-coordinate for a center of an i^(th) scanningspot;

[0018] (x_(0i)−s/2)≦x_(j)≦(x_(0i)+s/2);

[0019] y_(i)(x_(j)) is a depth of the i^(th) basis function for thegenerally straight treatment pattern; and

[0020] θ_(i)(x_(j)) is a coverage angle of the i^(th) basis function forthe generally circular treatment pattern.

[0021] In a specific embodiment, fitting the lens function with thebasis function comprises solving the following equation for coefficientsa_(i) representing treatment depth for the i^(th) scanning spot:${f\left( x_{j} \right)} = {\sum\limits_{i = 1}^{M}{a_{i}{X_{i}\left( x_{j} \right)}}}$

[0022] where

[0023] f(x_(j)) is the lens function; and

[0024] i=1, . . .l , M.

[0025] Fitting the lens function and the basis function may includespecifying a deviation for each of the N discrete evaluation points. Themethod may include refitting the lens function with the basis functionby varying the deviations to iterate for an acceptable fit or a bestfit. The method may further include refitting the lens function with thebasis function by varying the number of scanning spots M to iterate fora best fit. The method may also include refitting the lens function withthe basis function by varying the size of the scanning spot and/or thenumber of scanning pulses at a scanning spot location to iterate for abest fit. A merit function may be defined to represent an error of fitbetween the lens function and the basis function. The method includesminimizing the merit function to achieve the best fit. In someembodiments, the merit function and the total number of scanning spotsin the treatment table are both minimized to achieve the best fit withthe least number of scanning spots.

[0026] Another aspect of the invention is directed to a method ofgenerating scanning spot locations for ablating tissue using a scanninglaser beam for generating scanning spots over a treatment region largerin area than the scanning spots. The method comprises providing a lensfunction representing a desired lens profile for ablating the tissue byscanning spots of the laser beam on the tissue. A basis functionrepresents a treatment profile produced by the overlapping scanningspots along a treatment path. The basis function represents a sectionoriented across the treatment path. The lens function is fitted with thebasis function to obtain the scanning spot locations for the overlappingscanning spots of the laser beam. In specific embodiments, the treatmentprofile is symmetrical with respect to an axis of symmetry.

[0027] In accordance with another aspect of the invention, a method forfitting a three-dimensional target profile comprises providing atwo-dimensional basis function including overlapping portions torepresent a three-dimensional profile which has symmetry with respect toa two-dimensional section extending along a treatment pattern. Thethree-dimensional target profile is fitted with the two-dimensionalbasis function to obtain a distribution of the overlapping portions.

[0028] In accordance with yet another aspect of the present invention, asystem for ablating tissue comprises a laser for generating a laserbeam, and a delivery device for delivering the laser beam to a tissue. Acontroller is configured to control the laser and the delivery device. Amemory is coupled to the controller, and comprises a computer-readablemedium having a computer-readable program embodied therein for directingoperation of the system. The computer-readable program includes a firstset of instructions for generating a treatment table for ablating thetissue over a treatment region larger in area than the spot size of thelaser beam to achieve a desired lens profile for ablating the tissue, asecond set of instructions for controlling the laser to generate thelaser beam, and a third set of instructions for controlling the deliverydevice to deliver the laser beam to the tissue at the scanning spotlocations.

[0029] In specific embodiments, the first set of instructions includes afirst subset of instructions for providing a target functionrepresenting the desired lens profile for ablating the tissue byscanning spots of the laser beam on the tissue, a second subset ofinstructions for providing a basis function representing a treatmentprofile produced by the overlapping scanning spots in a treatmentpattern, and a third subset of instructions for fitting the targetfunction with the basis function to obtain the treatment table includingscanning spot locations and characteristics for the overlapping scanningspots of the laser beam.

BRIEF DESCRIPTION OF THE DRAWINGS

[0030]FIG. 1 is a top projection of a set of overlapping square pulsesillustrating a straight pattern;

[0031]FIG. 2 is a cross sectional projection of the overlapping squarepulses of FIG. 1 along the central meridian illustrating thecross-sectional treatment profile;

[0032]FIG. 3 shows a top projection of a treatment profile of a roundspot trench;

[0033]FIG. 4a is a schematic view of a circle rotated about the z-axis;

[0034]FIG. 4b shows a cross-sectional projection of an arbitrary shape;

[0035]FIG. 4c is a top projection of a set of overlapping circularpulses illustrating a ring pattern;

[0036]FIG. 5 shows the asymmetric treatment profiles of the overlappingcircular pulses in the ring pattern of FIG. 6 taken along a radialsection of the ring pattern;

[0037]FIG. 6 is a one-dimensional, simplified representation of a targetprofile fitting equation;

[0038]FIG. 7 is a flow diagram illustrating a simulated annealing methodfor solving a target profile fitting equation;

[0039]FIG. 8A is a block diagram of the laser system illustrating anembodiment of the present invention;

[0040]FIG. 8B is a block diagram of the control structure of a computerprogram according to a specific embodiment;

[0041]FIG. 9 shows the result of fitting a myopic lens equation using 27equally spaced basis functions;

[0042]FIG. 10 shows the result of an extended zone fit of a myopic lensequation with 47 equally spaced basis functions;

[0043]FIG. 11 is an interferometer image of a −4 diopter myopic lensablation profile using the simulation output from FIG. 10;

[0044]FIG. 12 shows ablation profile changes due to spot sizevariations;

[0045]FIG. 13 is an interferometer image of a hyperopic lens ablationprofile;

[0046]FIGS. 14A and 14B show extended hyperopic lens fit with twodifferent spot sizes;

[0047]FIG. 15 is an interferometer image of a −4D hyperopic cylinderablation profile;

[0048]FIGS. 16a-16 c are graphs illustrating the results of thesimulated annealing process at different stages of the iteration forfitting an elliptical profile;

[0049]FIGS. 17a and 17 b show contour plots of the elliptical targetprofile and the calculated profile; and

[0050]FIG. 18 shows a −4 diopter elliptical ablation profile based onthe solution obtained from the simulated annealing process of FIGS.16a-16 c.

DESCRIPTION OF THE SPECIFIC EMBODIMENTS

[0051] I. Generation of Scanning Locations and Depths

[0052] A. Treatment Space

[0053] Treatment space or profile is the treated form (or profile)created by a given spatial pulse sequence taking into account itsoverlap characteristics.

[0054] 1. Rectangular Spot Shape in Straight Pattern

[0055] As an introductory example, a simple square spot with a treatmentintensity of approximately 0.25 μm per pulse is used. The overlap is setto 33%. A 1 mm square spot centered at position (in mm) [0, 0] will havethe next spot at [0, 0.33]. The area the two pulses have in common is0.67 mm². The offset (step size) from the first pulse is 0.33 mm.Continuing this sequence for a total of N pulses will produce a treatedpath of width 1 mm, and length 0.33*N mm. The overlap is 33% so the lastpulse to coincide with the first pulse is the third pulse. The fourthpulse's edge is immediately adjacent to the first pulse. Assuming thatthe positioning is perfect, they will share a common edge. Looking infrom a meridian cross section, the slope is (depth per pulse)/0.33 mm.All positions after the third pulse will have a depth of (3*depth perpulse). The slope also occurs at the trailing end. The cut profile fromthe ends has the box shape, with a depth of (3 *depth per pulse) and awidth of 1 mm. The slope is:${Slope} = {\frac{{depth}\quad {per}\quad {pulse}}{{step}\quad {size}}.}$

[0056] The treatment profile of FIGS. 1 and 2 is generated by scanningoverlapping square spots in a straight treatment pattern. As shown inFIG. 1, the shape of the trench is rectangular in shape. The stacking ofpulses has the first interval with a single pulse depth, the secondinterval with a double pulse depth, and the third interval with a triplepulse depth, as seen in FIG. 2. The depth of the trench remains the sameafter the third pulse because the overlap is ⅓ of the spot, i.e., thefourth pulse does not overlap the first.

[0057] The treatment profile generated by overlapping rectangular spotsis relatively simplistic and easy to represent mathematically due to therectangular spot shape as well as the straight-line pattern. Thecomplexity of the mathematical representation increases for round spots,particularly if the spots overlap in a nonlinear pattern.

[0058] 2. Round Spot Shape in Straight Pattern

[0059]FIG. 3 shows a round spot having a diameter or spot size of s. xis a spatial coordinate from the origin to a point on the trenchproduced by the round spot (which may be discrete). x₀ is the positionof a basis function or the center of the ablation trench produced by theround spot, and θ is the coverage angle. A basis function is amathematical representation of the treatment profile produced, forinstance, by overlapping scanning spots in a specified treatmentpattern. Here, each round spot is assumed to produce an ablation with auniform ablation depth. Adjustments may be made for nonuniform energydistribution profiles, “central islands,” and other local ablation depthvariations, and the like.

[0060] Overlapping a round spot produces similar characteristics asoverlapping a rectangular spot. Repeating the same exercise foroverlapping the round spot in a straight pattern as done with the squarespot, the side profile will have the same look as the square spot. Theleading and trailing slopes will have the same slopes as those of thesquare spot. The profiles along the ends are different. It issemicircular in shape, or more precisely, semi-elliptical in shape, witha maximum depth of${Depth}_{\max} = {{\frac{{spot}\quad {size}}{{step}\quad {size}} \cdot {depth}}\quad {per}\quad {{pulse}.}}$

[0061] The slope is the same as described in the square spot. The roundspot has a different edge pattern. Without a better descriptive phraseto describe it, it has a “cloud” effect to its edge.

[0062] The basis function for a round spot in a straight pattern is:

y _(i)(x)={square root}{square root over ((spot size/2)²−(x−x _(0i))²)}

with

(x _(0i)−spot size/2)<x<(x _(0i)+spot size/2);

[0063] where

[0064] i is the running index of the basis functions;

[0065] x_(0i) is the position of the i^(th) basis function (i.e., thecenter of the i^(th) spot) measured from the origin;

[0066] y_(i)(x) is the trench depth of the i^(th) basis function; and

[0067] x is the spatial coordinate from the origin (e.g., optical axisof a cornea) to a point on the trench.

[0068] Fitting the function in a computer algorithm requires thecoordinate system to be discrete. Defining j as the running index forthe x spatial dimension, the trench depth is the same as above but withx replaced by x_(j),

y _(i)(x _(j))={square root}{square root over ((spot size/2)²−(x _(j) −x_(0i))²)}

[0069] where j is the running index for the sampling or evaluationpoints.

[0070] 3. Round Spot Shape in Ring Pattern

[0071] Generation of a ring treatment pattern makes use of axisymmetry.To illustrate how symmetry affects the overlap, consider a circular diskof radius r with its center b′ away from the origin as depicted in FIG.4a. Recall from differential geometry that its normal vector (right-handrule) gives the “direction of the surface”. The disk's direction vectorat anytime is in the x-y plane thus implying that the surface of thedisk itself lies in the y-z plane. Revolving this disk about the originin the direction of its normal vector (z-axis) will generate a bagelshaped volume. The bagel will have a radius of b′ and a thickness of 2r.This assumes that b′>r.

[0072] Consider again, the same disk but only the top half (z>0). Thebottom is cut off at the x-y plane. Rotating this half disk about thesurface along its normal vector will generate a semi “flat bottom slicedbagel” ring with a diameter of 2b′ and thickness of r. Integrating thevolume as the half disk moves about the origin will produce the volumeof the flat bottom bagel. This type of “differential volume” will allowus to produce a “bowl” shape (flat on one side, curved on the other)similar to the bottom half of a sliced bagel.

[0073] One method for calculating volume is a 3-D version of Simpson'srule. FIG. 4b shows a cross section of a set of circular rings at a, x₁,. . . , x_((n−1)) with a constant thickness of Δx. The top curve is thedesired target shape representing, for example, a lens volume. Thecolumns are cross sections of “thick” rings, spaced x_(j) apart. Summinga particular set of rings will produce the desired volume. Note thatthis is mathematically easy to do because the basis functions areorthogonal. The rings following Simpson's rule will need to be very thinto get a good approximation to the curve (i.e. the error is on the orderof θ²(x)).

[0074] Generalizing these approaches, symmetric objects can be generatedin a “space” minus one dimension. For example, a symmetric 3-D objectcan be generated in 2-D space if symmetry is utilized. Generating thesebasis “rings” or treatment volumes from differently sized circular spotsis the basis for one aspect of the present invention.

[0075] The thin lens volume is generated using flat bottom bagel shapedrings. Finding the different combinations of ring radii (b′) andthickness (r) to generate the lens is done via, for example, a linearleast square fitting routine.

[0076] Scanning a continuous circular spot in a large ring pattern willproduce a result that resembles a bagel with a half-circular crosssection (like a sliced bagel). As the ring's position gets closer to theorigin, the ring's cross section becomes very noncircular because theablation spot overlap characteristics change as the ring's radiuschanges.

[0077] The scan pattern for a round spot in a ring pattern will looksomewhat like a torus but with a radial dependence resulting in a bulgetoward the origin. To understand this radial dependence, consider a spotmoving along a ring as shown in FIG. 4c. The outer portion of the spotwill have an angularly smaller coverage than that of the inner portion.That is, the overlap for the outer portion of the spot is less than theinner portion. The functional expression for the arc length for thisspot as a function of the radial distance is an arccosine:${\theta_{i}(x)} = {\cos^{- 1}\left( \frac{x^{2} + x_{0i}^{2} - \left( {s/2} \right)^{2}}{2 \cdot x_{0i} \cdot x} \right)}$

[0078] or in terms of spot size or diameter,${\theta_{i}(x)} = {\cos^{- 1}\left( \frac{x^{2} + x_{0i}^{2} - \left( {{spot}\quad {{size}/2}} \right)^{2}}{2 \cdot x_{0i} \cdot x} \right)}$

[0079] where

[0080] i is the running index of the basis functions;

[0081] x is the spatial coordinate from the origin (e.g., optical axisof a cornea) to a point on the trench; and

[0082] x_(0i) is the position of the i^(th) basis function (i.e., thecenter of the i^(th) spot) measured from the origin; and

[0083] θ_(i)(x) is the coverage angle of the i^(th) basis function.

[0084]FIG. 3 shows that the i^(th) basis function is defined fromx_(0i)−s/2<x<x_(0i)+s/2 and zero otherwise. x is continuous but actualimplementation in an algorithm to solve for the scanning coordinates orlocations will have discrete x.

[0085] A discrete coordinate system facilitates fitting the function ina computer algorithm. Defining j as the running index for the x spatialdimension, the angle of coverage is the same as above but with xreplaced by x_(j),${\theta_{i}\left( x_{j} \right)} = {\cos^{- 1}\left( \frac{x_{j}^{2} + x_{0i}^{2} - \left( {{spot}\quad {{size}/2}} \right)^{2}}{2 \cdot x_{0i} \cdot x_{j}} \right)}$

[0086] where j is the running index for the sampling or evaluationpoints.

[0087] A feature of these treatment “rings” is the asymmetric treatmentspace profile with respect to the radial direction from the origin asshown in FIG. 5. The profiles get more asymmetric as the pulse ringsapproach the origin, as illustrated by the increased asymmetry ofprofile A over profile B. This feature allows the fitting of a centrallymaximum curved surface. The coefficient of each basis function will givethe treatment depth.

[0088] 4. Treatment Matrix

[0089] The above examples employ specific treatment profiles that can beexpressed in analytical forms for achieving relatively simple types oftarget ablation profiles. Other treatment profiles may be used. Ageneral treatment matrix includes data of treatment profiles produced byoverlapping scanning spots of any desired shapes and sizes.

[0090] At one time it was thought that the laser-tissue interactionresulted in a flat ablation profile, so that each pulse could be thoughtof as removing a disc of material, with a diameter equal to the laseriris width and a uniform depth. Empirical evidence, however, has shownthat the laser-tissue interaction may result in a rather complicatedablation profile with a complex shape and characteristics that can varysignificantly with the laser iris diameter. In general, the ablationprofiles for different laser spot sizes can be digitized and arranged asa treatment matrix of ablation data to be used for numericallydetermining a treatment plan or table (i.e., ablation amplitudes for atreatment pattern of scanning spots) to achieve a target ablationprofile. There will be a different set of data for each treatment mediumincluding, for example, polymethylmethacrylate (PMMA), calibrationplastic, human cornea, porcine cornea, or the like.

[0091] B. Target Profile Fitting

[0092] Fitting the basis function to the target function representingthe target ablation profile produces the scanning spot locations andablation depths. In specific embodiments, the target function isrepresented by the lens equation as described below.

[0093] The target profile fitting involves solving for a set ofcoefficients for the following target profile fitting matrix equation:${f(x)} = {\sum\limits_{i = 1}^{M}{a_{i}{X_{i}(x)}}}$

 F=AX

[0094] with the solution as an inversion

A=X ⁻¹F

[0095] where

[0096] X_(i)(x) is the i^(th) basis function;

[0097] f(x) is a particular target profile, which may be expressed as alens equation or represented by a target matrix of data; and

[0098] a_(i) is the amplitude of the i^(th) basis function.

[0099] F is the lens equation or target matrix, X are the basisfunctions, and A is the treatment plan or table containing theamplitudes for a particular treatment pattern (i.e., the scanning spotlocations and characteristics such as size, shape, and depth). Thefitting program calculates A such that the error (χ²'s) between F and Xis a minimum.

[0100] Some of the examples discussed above employ analytical basisfunctions. For instance, the basis functions X are y_(i)(x_(j)) for thestraight pattern and θ_(i)(x_(j)) for the ring pattern. The straightpattern is applicable in myopic and hyperopic cylinders, while the ringpattern is applicable in myopia and hyperopia. For phototherapeutickeratectomy (PTK), the treatment profile is essentially a centrally flathyperopia ablation since PTK is a flat ablation. Thus, the basisfunction θ_(i)(x_(j)) can be used for PTK (while the lens equation forPTK is a constant value). Of course, the basis function may berepresented by a general treatment matrix instead of the analyticalfunctions. Furthermore, the target profilef(x) may be expressed as alens equation which has been developed for certain cases such as myopiaand hyperopia. Other target ablation profiles (e.g., an ellipticalprofile, a toric profile, a wavefront profile, or an arbitrary profile)may be expressed more generally as a target matrix.

[0101]FIG. 6 schematically illustrates in one dimensional simplifiedform the equation XA=F, where X is a basis function expressed in ageneral treatment matrix representing ablation profiles that may haverelatively complex shapes, F is the target profile schematically shownas a lens, and A is the treatment plan or table to be solved. In thisexample, the treatment table is a pulse instruction vector (PIV)representing the laser instruction count for each basis function. In thespecific example shown, the treatment matrix represents ablationprofiles for five laser pulses (i.e., five particular laser pulsediameters or shapes at particular locations). The elements of thetreatment table or PIV correspond to the numbers of pulses that can becombined to build an entire ablation having a shape that matches thetarget profile as close as possible. In that case, the elements of thetreatment table or PIV are positive integers. The treatment table or PIVis used to operate the laser to achieve ablation that approximates thetarget profile. Note that a two-dimensional version of the equation(instead of the one-dimensional illustration shown in FIG. 6) willproduce a three-dimensional model of the ablation profile.

[0102] 1. Least Square Fitting Approach Using Analytical Lens Equations

[0103] The difficulty in performing the inversion of the target profilefitting equation to obtain A is that X⁻¹ is not orthogonal. One usefulmethod for solving this equation is to use a least square fittingroutine. An example of a suitable method for fitting the basis functionis found in Numerical Recipes in C, 2^(nd) Ed., W. H. Press, S. A.Teukolsky, W. T. Vetterling, B. P. Flannery, Cambridge Univ. Press, pp.671-675 (1992). The general model representation is${f(x)} = {\sum\limits_{i = 1}^{M}{a_{i}{X_{i}(x)}}}$

[0104] where X₁, . . . , X_(M) are arbitrary functions of X, calledbasis functions, and M is the number of basis functions. The X_(i)'sneed not be linear. The linearity is in respect to the linear dependenceon the parameters a_(i). The basis functions are the y_(i)'s for thestraight treatment pattern or the θ_(i)'s for the ring treatmentpattern.

[0105] The merit function (χ²) is defined as$\chi^{2} = \left\lbrack \frac{{f\left( x_{j} \right)} - {\sum\limits_{i = 1}^{M}{a_{i}{X_{i}\left( x_{j} \right)}}}}{\sigma_{i}} \right\rbrack^{2}$

[0106] f(x) represents the target profile to be fitted with the basisfunctions. In specific embodiments, f(x) is the lens equation evaluatedat x_(j), that is, f(x_(j)). The σ_(i)'s are the uncertainties orstandard deviations of the X_(i)'s. These are typically set to one butcan be altered if a more precise fit for certain parts of the curve isdesired (e.g., the central half can be set to 0.25 and the outer halfcan be set to 1). For instance, decreasing the value of σ_(i) for thecentral area and increasing it for the outlying area will force the fitto be more precise in the central region. The reduced χ² is$\chi_{\upsilon}^{2} = \frac{\chi^{2}}{\upsilon}$

[0107] where υ is the number of degrees of freedom.

[0108] To minimize χ², the derivative of χ² with respect to parametersa_(k) for all M's is taken and set to zero.$0 = {\sum\limits_{i = 1}^{M}{{\frac{1}{\sigma_{i}^{2}}\left\lbrack {{f\left( x_{j} \right)} - {\sum\limits_{i = 1}^{M}{a_{i}{X_{i}\left( x_{j} \right)}}}} \right\rbrack}{X_{k}\left( x_{j} \right)}}}$

[0109] for i=1, . . . , M. These are the normal equations of the leastsquare problem. The two commonly used methods of obtaining the vector ais by LU decomposition and by Gauss-Jordon elimination. The Gauss-Jordonmethod will provide the variances of a_(i) as part of the solution.

[0110] χ_(υ) ² gives the relative “goodness” of fit. The better the fit,the closer χ_(υ) ² approaches one (χ_(υ) ²→+1). If χ_(υ) ² is less thanone, it means the original estimates of σ_(i) ² 's are too large.

[0111] The input parameters used by the fitting routine are:

[0112] (1) The data generated by the lens equation, f(x_(j)).

[0113] (2) The standard deviation σ_(i) of each data point. Since thedata is generated from the lens equation, the standard deviation isarbitrarily set to one. This arbitrariness needs to be addressed becausethe standard deviation affects χ². Consistency is desired because oncethe values for the standard deviation is determined, it should not bechanged if χ² is to be used to determine the relative goodness of fit. Awide variety of fit evaluation methods could alternatively be usedwithin the scope of the invention.

[0114] (3) The Basis Functions in Analytical Form.

[0115] 2. Simulated Annealing Approach

[0116] Another computational technique that may be used for solving thetarget profile fitting equation is simulated annealing. Such a techniqueis described, for instance, in Numerical Recipes in C, 2^(nd) Ed., W. H.Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, CambridgeUniv. Press, pp. 444-451 (1992). Simulated annealing is an evolutionaryoptimization method using gradual, randomly induced changes to find thebest solution to mathematically intractable problems.

[0117] Simulated annealing is a method used for minimizing (ormaximizing) the parameters of a function. It is particularly suited toproblems with very large, poorly behaved function spaces. Simulatedannealing can be applied in the same way regardless of how manydimensions are present in the search space. It can be used to optimizeany conditions that can be expressed numerically, and it does notrequire a derivative. It will not be fooled by local minima in thesearch space.

[0118] Simulated annealing is inspired by the observation that if athermodynamic system is cooled slowly, it will end up very near itsminimum energy state. For example, a piece of hot metal contains nmolecules in rapidly changing orientations. If the metal is quenched incold water the thermal energy rapidly leaves the metal and the moleculesremain—“frozen forever”—in their last orientation. However, if the metalis allowed to cool gradually the molecules will continue changingorientation, but gradually favoring lower and lower energy states.Eventually all the molecules will be “frozen” into their lowest energystate. (In metals, the molecules align, forming crystals.) This iscalled annealing.

[0119] Instead of a piece of metal, the treatment table or pulseinstruction vector (PIV) will be annealed. The PIV may be expressed as a1 by N array of integers, each of which describes the number of pulsesat a particular location and given pulse size, shape, and/or depth. ThePIV may be expressed as follows: $\begin{bmatrix}{\# \quad {of}\quad 1\quad {mm}\quad {pulses}\quad {at}\quad {position}\quad 1} \\{\# \quad {of}\quad 1\quad {mm}\quad {pulses}\quad {at}\quad {position}\quad 2} \\{\# \quad {of}\quad 1\quad {mm}\quad {pulses}\quad {at}\quad {position}\quad 3} \\\cdots \\{{\# \quad {of}\quad 6\quad {mm}\quad {pulses}\quad {at}\quad {position}\quad M} - 1} \\{\# \quad {of}\quad 6\quad {mm}\quad {pulses}\quad {at}\quad {position}\quad M}\end{bmatrix}\quad$

[0120] This vector can be very large. Some examples have upwards of 3000elements.

[0121] To “simulate” annealing, two physical concepts are applied tosolving the target profile fitting equation. The first is the concept ofenergy state to be minimized, which is analogous to the error or meritfunction (χ²) defined as$\chi^{2} = {\sum\limits_{i = 1}^{M}\left\lbrack {{f\left( x_{j} \right)} - {a_{i}X_{i}}} \right\rbrack^{2}}$

[0122] where f(x) represents the target profile to be fitted with thebasis functions X_(i). In specific embodiments, f(x) is the targetprofile evaluated at x_(j), that is, f(x_(j)). The numerical value (χ²)to be minimized is the sum of the squares of the difference at eachpoint between the target shape and the “test shape” or “working shape”during the iteration. This example employs a least square fit forminimizing the merit function.

[0123] The merit function may be modified to find a fit with moredesirable optical qualities. For instance, a weighted merit function maybe$\chi^{2} = {\sum\limits_{i = 1}^{M}\left\lbrack \frac{{f\left( x_{j} \right)} - {a_{i}X_{i}}}{\sigma_{i}} \right\rbrack^{2}}$

[0124] The σ_(i)'s are the uncertainties or standard deviations of theX_(i)'s. These are typically set to one but can be altered if a moreprecise fit for certain parts of the curve is desired (e.g., the centralhalf can be set to 0.25 and the outer half can be set to 1). Forinstance, decreasing the value of σ_(i) for the central area andincreasing it for the outlying area will force the fit to be moreprecise in the central region.

[0125] The second physical concept that is applied in the fittingapproach is that of temperature, which is analogous to a measure of therandom change applied to the treatment table or PIV (i.e., change in thenumber of laser pulses) before each successive evaluation. In this case,each “degree” of “temperature” change will cause one random fluctuationin the PIV. For each fluctuation one element of the PIV will be randomlyselected. That element will then be either incremented or decremented,again based on a random number. In the case that the element is alreadyzero and a decrement is chosen, nothing will happen to that element,since only positive integers are acceptable. Evaluation of the “energystate” (i.e., merit function) takes place during the “temperature”change iteration for the PIV until a sufficiently low value for themerit function is achieved.

[0126] An example of implementing the simulated annealing procedure oralgorithm is illustrated in FIG. 7. After the master application callsthe algorithm (step 20) and the simulated annealing subroutine isinitiated (step 22), the subroutine randomly modifies the treatmenttable or pulse interaction vector (PIV) according to the currenttemperature (i.e., size of the incremental change to an element of thePIV) to create the next vector PIV′ (step 25). An evaluation takes placein step 26 to determine whether PIV′ is better than PIV in reducing themerit function or “energy state.” If PIV′ is better than PIV, then thevector is updated by setting it to PIV′ (step 28), and the iterationcontinues by following the loop back to step 24. If PIV′ is not betterthan PIV, then the subroutine determines whether the vector has remainedconstant for a preset number (n) of evaluations (step 30). If not, thePIV is not updated and the iteration continues by returning to step 24.If n number of evaluations have taken place with no improvement of thePIV, then the subroutine determines in step 32 whether the temperature(i.e., the size of the incremental change) has reached a preset minimumvalue. If not, the temperature is lowered by one degree in step 34(i.e., the size of the incremental change to an element of the PIV islowered by a preset amount) and the iteration continues by returning tostep 22. If the minimum temperature has been reached, then the annealingis complete (36) and the simulated annealing procedure is terminated.Therefore, the annealing process constitutes repeated evaluation of thePIV, at gradually decreasing temperatures. For a particular temperaturethe iteration cycle will repeat until a certain n number of evaluationshave failed to show improvement. Once the solution meets somepre-determined “goodness of fit criteria” as represented by the meritfunction, the final solution is returned.

[0127] For problems with radial symmetry, such as spheres or roundPTK's, the problem is solved in only one dimension. For problems withbilateral symmetry, such as ellipses or cylinders, the problem becomestwo dimensional. The method only needs to be used to solve one quarterof the surface and then the solution is applied over all four quadrants.This saves considerable computation. For problems with no symmetry, suchas a wavefront map, the method solves for the entire surface.

[0128] Solving the target profile fitting equation numerically for ageneral two-dimensional surface, particularly without the benefit ofsymmetry, can be very difficult and can require a great deal ofcomputational time. The simulated annealing approach can generallyproduce the best solution in a reasonable amount of time, preferably inreal time.

[0129] The input parameters used by the simulated annealing routine are:

[0130] (1) The target profile data, which can represent any arbitraryshape. The target shape can be arbitrary, allowing for unlimitedablation possibilities, including the map from a wavefront type deviceor any other topography system. As long as the target profile can bedescribed mathematically or with an array of data, the algorithm can useit. Radial or bilateral symmetry is not required, although the presenceof symmetry can make the computation faster.

[0131] (2) The treatment matrix (i.e., basis function) which canrepresent any ablation sizes and shapes. The laser ablation shape doesnot have to be flat or smooth, or have specific shapes. As long as itcan be measured and described mathematically or with an array of data,the simulated annealing process can be used. A unique description isprovided for every unique pulse shape or size to be used. A fixedspot-size laser would have only one description, while a variablespot-laser could have as many as desired.

[0132] (3) The standard deviation σ_(i) of each data point if a weightedmerit function is used. It is noted that any merit function can be used.This allows the user to decide what aspect of the solution is mostimportant. For example, extra “weighting” may be provided in the centralportion of the lens, so that the algorithm will recognize that thecenter portion of the lens is the most important part and will sacrificeaccuracy around the edges to get the best possible fit in the centerportion. Moreover, the merit function is not limited to just the shapeof the surface. For instance, the merit function can be linked to thenumber of pulses, so that the algorithm will favor solutions with fewerpulses. This helps to minimize the treatment time. In one example, themerit function is expressed as follows:$\chi^{2} = {\left\{ {\sum\limits_{i = 1}^{M}\left\lbrack \frac{{f\left( x_{j} \right)} - {a_{i}X_{i}}}{\sigma_{i}} \right\rbrack^{2}} \right\} \times \left( {\left( N_{{number\_ of}{\_ pulses}} \right)^{2} + K_{offset}} \right)}$

[0133] In this example, the merit function is just the surface fitfunction multiplied by the number of pulses squared plus a constantoffset. The offset lets the user tune the relative importance of thesurface fit versus the pulse count. If the offset is large compared tothe number of pulses, then the merit function will be only mildlyinfluenced by pulse count. The solution will tend to reduce the numberof pulses, but only slightly. If the offset is comparable or smallcompared to the number of pulses then the merit function will besignificantly influenced. The solution will compromise some goodness offit in order to use fewer pulses. It is possible to have solutions withdiffering pulses because a variety of pulse sizes are used (e.g., asingle five millimeter pulses can remove as much tissue as twenty fivepulses of one millimeter diameter).

[0134] Other modification of the merit functions would includemaximizing smoothness, (i.e. minimizing the first derivative of thesurface error), meeting hardware constraints of the Laser deliverysystem (e.g., the Laser can only move so far in a given time), meetingbiological constraints on the eye (e.g., pulses might be distributedover space and time to avoid excess thermal build up.)

[0135] C. Software Implementation Steps

[0136] 1. Least Square Fitting Approach Using Analytical Lens Equations

[0137] The following steps are provided to illustrate an example ofsoftware implementation of the fitting approach of generating scanninglocations and ablation depths using analytical lens equations.

[0138] Step 1. Evaluate Lens Equation. The lens equation is evaluatedfor N points (typically equally spaced) in the region 0≦x<L where L isthe treatment zone radius for a ring pattern and is the treatment zonelength for a straight pattern. N is set to about 300 in specificembodiments. If it is desired to emphasize the central portion of thelens equation (e.g., for myopia and myopic cylinder), a “shift” issubtracted from the treatment zone length.

[0139] For myopia and myopic cylinder, the lens equation is:${f\left( x_{j} \right)} = {\sqrt{R_{1}^{2} - x_{j}^{2}} - \sqrt{\left( \frac{R_{1}\left( {n - 1} \right)}{n - 1 + {R_{1}D}} \right) - x_{j}^{2}} + C}$

[0140] where

[0141] 0≦x_(j)≦(L−shift); and

[0142] j=0, 1, . . . , N−1.

[0143] The constant C is defined as$C = {\sqrt{R_{1}^{2} - {s^{2}/4}} + \sqrt{\left( \frac{R_{1}\left( {n - 1} \right)}{n - 1 + {R_{1}D}} \right) - \frac{s^{2}}{4}}}$

[0144] where

[0145] x_(j) is the reference x-coordinate measured from the opticalaxis of a cornea of the j^(th) evaluation point for the center of theablation or scanning spot;

[0146] s is the diameter of the ablation (i.e., spot size);

[0147] R₁ is the anterior radius of curvature of the cornea in meters;

[0148] R₂ is the final anterior radius of curvature of the cornea inmeters;

[0149] n=1.377 is the index of refraction of the cornea;

[0150] D is the lens power of the ablation in diopters;

[0151] L is the treatment zone length representing the length across agenerally uniform section which is oriented normal across a generallystraight treatment pattern or is oriented radially across a generallycircular treatment pattern; and

[0152] shift is the amount of emphasis shift (shift typically is about0-0.2, and can be set to zero).

[0153] Note that the above lens equation applies for the full range ofangles in myopia to provide an axisymmetric profile, but applies only ata particular angle in myopic cylinder to provide a constant profilealong a straight pattern.

[0154] For hyperopia and hyperopic cylinder, the lens equation is:${f^{\prime}\left( x_{j} \right)} = {R_{1} - \frac{R_{1}\left( {n - 1} \right)}{n - 1 + {R_{1}D}} - \sqrt{R_{1}^{2} - x_{j}^{2}} + \sqrt{\left( \frac{R_{1}\left( {n - 1} \right)}{n - 1 + {R_{1}D}} \right) - x_{j}^{2}}}$

[0155] where

[0156] 0≦x_(j)≦(L−shift); and

[0157] j=0, 1, . . . , N−1.

[0158] Note also that this lens equation applies for all angles inhyperopia and at a particular angle in hyperopic cylinder. A moredetailed description of the lens equation is provided in Charles R.Munnerlyn et al., “Photorefractive Keratectomy: A Technique for LaserRefractive Surgery,” J. Cataract Refract. Surg. (January 1988), which isincorporated herein by reference in its entirety.

[0159] For PTK, the lens equation is:

f′(x _(j))=d

[0160] where d is a constant depth.

[0161] Step 2. Define basis rings at different radial positions for aring pattern or define basis strips at different lateral positions for astraight pattern. A set of radial positions for basis rings or a set oflateral positions for basis strips are defined at x_(0i). In specificembodiments, the number of rings/strips, M, may vary from about 17 toabout 57. The positions of the ring/strip's radial center are equallyspaced points in the region 0≦x≦(L−spot size/2+extended zone). Theextended zone is for the taper or “feathering” of the outercircumference or edge. The extended zone typically ranges from about0.1-0.5 mm, and may be fixed at about 0.2 mm. The extended zone alsohelps in the fitting routine. M is the number of basis functions. Thenumber should be large enough to give a smooth fit and not too large tocreate quantized errors. Typically, M is about 7-97. In specificembodiments, Mstarts at about 17 for -1.5D or less, becomes about 27between −1.5 and -2.5 D, and reaches about 47 for greater than -2.5 D.The i^(th) basis function is

x _(0i)=[(treatment zone length−spot size+extended zone)/M]*i

[0162] wherein i=1, 2, . . . , M.

[0163] Step 3. Define the positions for lens equation evaluations. N isthe number of evaluation points. The i^(th) position for lens equationevaluation is

x _(j)=[treatment zone length−shift)/N]*j

[0164] where j=0, 1, . . . , N−1.

[0165] Step 4. Assign a deviation, σ_(j), for each x_(j). The deviationswill determine how tight the fit will be for each point. In a specificembodiment, the deviations are set to about 0.25 for the first 200points as 0.25 and about 1.0 for the final 100 points. The deviation canbe changed and optimized for different profiles and patterns.

[0166] Step 5. Apply least square fit. The function to fit the basisfunctions to the lens equation (i.e., fit X(x_(j)) to j(x_(j))) is${f\left( x_{j} \right)} = {\sum\limits_{i = 1}^{M}{a_{i}{X_{i}\left( x_{j} \right)}}}$

[0167] where, for a ring pattern (myopia, hyperopia, and PTK),${X_{i}\left( x_{j} \right)} = {{\theta_{i}\left( x_{j} \right)} = {\cos^{- 1}\left( \frac{x_{j}^{2} + x_{0i}^{2} - \left( {{spot}\quad {{size}/2}} \right)^{2}}{2 \cdot x_{0i} \cdot x_{j}} \right)}}$

[0168] and, for a straight pattern (myopic cylinder and hyperopiccylinder),

X _(i)(x _(j))=y _(i)(x _(j))={square root}{square root over ((spotsize/2)²−(x _(j) −x _(0i))²)}

[0169] Note that cos⁻¹ is valid only in the range (x_(0i)−spotsize/2)≦x_(j)≦(x_(0i)+spot size/2) for the ring pattern. The index j isthe running index for evaluation and fitting points. The index i is therunning index for the number of basis rings or basis strips.

[0170] The least square fit routine is an object that receives x_(j)coordinates, the lens equation evaluated at these coordinates (f(x_(j))values), the basis functions X_(i)(x_(j)), and returns a row (or column)matrix A representing the coefficients, and a χ². Once again, a varietyof fit or optimization routines might alternatively be employed.

[0171] Step 6. Evaluate closeness of fit. χ² is recalculated ifnecessary, that is, if any changes are made to the returned matrix A. χ²will give the goodness of fit. The difference between the lens equationand sum of the basis functions (standard deviation, etc.) can beevaluated to check for cusps. Cusps are small deviations with sharptransitions typically occurring at an integer multiple of the spot size.The cusps tend to be “stiff” equations and are difficult to eliminatelater. Cusps tend to be localized, occurring within 20 or so coordinatepoints. A moving average and standard deviation “window” of about 20points can be done to test for cusps. Note that the effects near theedge of the ablation will typically have a large standard deviation. Ifthe deviation due to cusps is too large, different fitting parameterscan be used to re-evaluate the closeness of the fit. The results shouldbe checked for negative coefficients and their magnitude. Another fitwith a different set of rings (usually a smaller number but may be alarger number) will need to be done if zeroing out the negativecoefficients causes a large χ². For small values, it suffices to set thenegative coefficients to zero and re-evaluate χ² and check the lensequation versus the solution.

[0172] Step 7. Iterate for best solution. If there are 2 cusps (criteriato be determined) or if χ² is too large (>100 for the initial settingsof 0.25 and 1.0), or if the solution is not close to the lens equation(to be determined), refitting is performed using different fitparameters:

[0173] (A) Begin with the number of basis functions, searching from17-57 at the initial spot size (2 mm);

[0174] (B) Next search with a different spot size (1.75 mm, 1.5 mm) andrepeat (A);

[0175] (C) Next search with extended zone and treatment size such thatthe sum of extended zone and treatment dimensions is a constant (i.e.,6.2 mm) in 0.05 mm increments, and repeat (A) and (B); and

[0176] (D) Next search with different σ_(j)'s. This can get rathercomplicated. In one calculation, the initial 200 σ_(j)'s were 0.25, andthe final 100 σ_(j)'s were 1.00. The initial 200 σ_(j)'s were changed to0.5 and (A), (B), and (C) were repeated. Then the 200 σ_(j)'s werechanged to 199 and (A), (B), and (C) were repeated, and so on.

[0177] Step 8. Determine actual pulse locations. The coefficients aireturned from the fit (matrix A) are used to calculate the number ofpulses at each radial position for a ring pattern or each lateralposition for a strip pattern (i.e., basis functions),

n _(k) =a _(k)/(depth per pulse)

[0178] With the present fitting routine, if the coefficient is less thanzero, it is typically set to zero. For myopia, the coefficients arerarely zero.

[0179] For a ring pattern, the pulses are spaced evenly along thecircumference, but with a random angle added to the first pulse.Dividing this number into 2π will give the angular spacing between eachpulse. The first pulse of each ring is desirably randomized to avoid apattern effect.

[0180] Step 9. Randomize table. All pulses are grouped into one andrandomized before treatment to avoid sequentially generating spatiallyadjacent treatment pulses.

[0181] Geometrical considerations can be addressed during the fittingprocess for a ring pattern (as in myopia or hyperopia) or after thefitting process for a straight pattern (as in myopic and hyperopiccylinders). For reference, the equation for the treatment space profilefor the ring pattern in the fitting step is repeated here:${a_{i}{\theta_{i}\left( x_{j} \right)}} = {a_{i}{\cos^{- 1}\left( \frac{x_{0i}^{2} + x_{j}^{2} - \left( {{spotsize}/2} \right)^{2}}{2 \cdot x_{0i} \cdot x_{j}} \right)}}$

[0182] The a_(i)'s are the basis function amplitudes returned by thefitting function (step 7). Each amplitude represents the height of a“ring”. The length of the ring is 2πx_(0i), giving the total number ofpulses as for the ring as follow:${{total}\quad {pulses}} = {\frac{2\quad \pi \quad x_{0i}}{{spot}\quad {size}} \cdot \frac{a_{i}}{{depth}\quad {per}\quad {pulse}}}$

[0183] Because the basis functions have radial dependencies, it iseasier to implement the normalization factor during the fit by using thenormalized basis functions in the fitting step (step 5) as follows:${{\overset{\_}{a}}_{i}{{\overset{\_}{\theta}}_{i}\left( x_{j} \right)}} = {{\frac{{\overset{\_}{a}}_{i}}{\pi}{\theta_{i}\left( x_{j} \right)}} = {\frac{{\overset{\_}{a}}_{i}}{\pi}{\cos^{- 1}\left( \frac{x_{0i}^{2} + x_{j}^{2} - \left( {{spot}\quad {{size}/2}} \right)^{2}}{{spot}\quad {{size} \cdot x_{j}}} \right)}}}$

[0184] The {right arrow over (a)}_(i)'s returned from these basisfunctions are normalized. The total number of pulses per ring is${{total}\quad {pulses}} = \frac{{\overset{\_}{a}}_{i}}{{depth}\quad {per}\quad {pulse}}$

[0185] The pulses are spread evenly along the 2π arc. According to aspecific embodiment of randomizing the pulses for treatment in step 9,random phases are added to the starting points (the first pulse perradii). The random phases will reduce the periodicity caused by the“finiteness” of the pulse number.

[0186] 2. Simulated Annealing Approach

[0187] The steps for implementing the method of generating scanninglocations and ablation depths utilizing the simulated annealing approachmay be analogous to those for the method that employs the least squarefitting approach using analytical lens equations.

[0188] Step 1. Provide target profile matrix. The target profile matrixcontains data that may represent any arbitrary target profile, or may beobtained by evaluating the lens equation at discrete points.

[0189] Step 2. Define the basis function in the form of treatmentmatrix. The treatment matrix contains data that represent ablationprofiles (i.e., sizes and shapes). The ablation profiles can bearbitrary as long as they can be measured.

[0190] Step 3. Define the positions for evaluations. The positions forevaluations may be selected, for instance, based on the target profile,the ablation profiles, or the like.

[0191] Step 4. Define the merit function. Any desirable merit functionmay be specified. An example is a weighted least square fit with astandard deviation σ_(i) specified for each data point.

[0192] Step 5. Apply simulated annealing. The simulated annealingprocess as illustrated in FIG. 7 and described above is applied to fitthe treatment matrix to the target profile matrix to solve for thetreatment plan or table containing the amplitudes for a treatmentpattern which is in general a two-dimensional pattern.

[0193] Step 6. Evaluate closeness of fit. The closeness of the fit ismonitored by evaluating the merit function which can be a least squarefit error function or any arbitrary function specified by the user. Theevaluation occurs in step 26 in the flow chart shown in FIG. 7.

[0194] Step 7. Iterate for best solution. As illustrated in FIG. 7, theiteration scheme involves random variation of the PIV and evaluation ofthe merit function (steps 24-30), and adjusting the “temperature”, i.e.,the size of incremental change to an element of the PIV (steps 32 and34).

[0195] Step 8. Determine actual pulse locations. The treatment plan ortable obtained by solving the target profile fitting equation includesamplitudes of the pulses which, in conjunction with the treatmentmatrix, provides instructions to the laser to direct pulses with thesizes and depths at the desired locations to achieve the best fit to thetarget profile.

[0196] Step 9. Randomize table. The pulses are randomized for treatmentto reduce the periodicity caused by the finiteness of the pulse number.

[0197] II. Exemplary Scanning System

[0198]FIG. 8A shows a block diagram of a scanning system 100 forimplementing the method of the present invention to perform cornealsculpting. The scanning system 100 includes a laser 102 for generating alaser beam and laser delivery system optics 104 for delivering the laserbeam to a target which is the cornea 106 of an eye. The laser deliverysystem optics 104 include various optical elements used to focus,modify, and direct the beam in a scanning mode to the cornea 106. Theenergy level of the laser beam is typically about 50-250 mJ/pulse.Commercially available systems having suitable laser and optics forscanning include, for example, the VISX TWENTY/TWENTY EXCIMER LASERSYSTEM, from VISX, Incorporated of Santa Clara, Calif.; Lasersight,Orlando, Fla.; Chiron Technolas, Germany; and Autonomous Technology,Orlando, Fla.

[0199] Optionally, the system 100 may include a calibration tool 108 fordetermining the characteristics of the laser beam which can be used asfeedback for calibrating the scanning system 100. In one embodiment, thecalibration tool 108 includes a reference-edge and a photodetector. Thecalibration tool 108 is placed at a target location for sensing a laserbeam directed onto the tool 108. The laser beam can be repeatedlyre-directed between the tool 108 and the patient's cornea 106. As such,after determining the size, shape, and/or position of the beam using thetool 108, the laser beam can be applied at a known location on thecornea. A repetitive measurement of intensity and shape characteristicsof the laser beam can be made, and a repetitive recalibration of thetargeting of the laser beam can be achieved to ensure precise positionalaccuracy when ablating the cornea 106. Using the sensed information fromthe calibration tool 108, an algorithm for calculating the locations andnumber of scanning spots (such as that according to the presentinvention) can be revised, thereby increasing the accuracy of thesculpting procedure. This calibration information can be used to adjustthe ablation algorithm immediately before and/or during each ablationprocedure in real time. A more thorough discussion of the calibrationprocedure is found in commonly assigned U.S. Patent Application entitled“Method and Apparatus for Determining Characteristics of a Laser BeamSpot,” Ser. No. 09/395,809 (Attorney Docket No. 18158-124), filed Sep.14, 1999, which is incorporated herein by reference in its entirety.

[0200] The scanning system 100 further includes a controller 110 forcontrolling operation of the laser 102 and the delivery system optics104, and all activities of the system 100. In one embodiment, thecontroller 110 includes a hard disk drive (memory 112), a floppy diskdrive, and a processor 114. The controller 110 executes system controlsoftware, which is a computer program stored in a computer-readablemedium such as the memory 112. The memory 112 is typically a hard diskdrive, but may also be other kinds of memory. The computer programincludes sets of instructions that dictate, for instance, spot ablationpattern, spot ablation depth, depth per pulse, the energy level, spotsize and shape, pulse rate, pulse location, and scanning pattern andsequence. It is understood that other computer programs stored on othermemory devices including, for example, a floppy disk or anotherappropriate drive, may also be used to operate the controller 110.

[0201] An interface 120 is provided between a user and the controller110 typically in the form of a display monitor for displayinginformation, and an input device such as a keyboard, a mouse, and/or alight pen to allow the user to communicate with the controller 110.

[0202] The scanning process can be implemented using a computer programproduct that is executed by the controller 110. The computer programcode may be written in any conventional computer readable programminglanguage. Suitable program code is entered into a single file, ormultiple files, using a conventional text editor, and stored or embodiedin a computer usable medium, such as a memory system of the computer. Ifthe entered code text is a high level language, the code is compiled,and the resultant compiler code is then linked with an object code ofprecompiled library routines. To execute the linked, compiled objectcode, the system user invokes the object code, causing the computersystem to load the code in memory. The processor 114 then reads andexecutes the code to perform the tasks identified in the program.

[0203]FIG. 8B shows an illustrative block diagram of the controlstructure of the system control software, computer program 130,according to a specific embodiment. Through the input device of theinterface 120, a user enters a set of process parameters into a processcontrol subroutine 132 in response to menus or screens displayed on themonitor. The process parameters include input data needed to operate thesystem 100 including, for example, the laser beam energy and spot size,and all the input necessary to determine the pulse locations and numbersaccording to the algorithm described above.

[0204] A scanning location determination subroutine 134 includes programcode for accepting the process parameters from the process controlsubroutine 132, and for determining the scanning pulse locations andnumbers according to the algorithm described above. The scanninglocation determination subroutine 134 produces laser control parametersand optics control parameters used in operating the laser 102 anddelivery system optics 104. A laser control subroutine 136 includesprogram code for accepting laser control parameters to control operationof the laser 102. An optics control subroutine 138 includes program codefor accepting optics control parameters to control operation of thesystem optics 104. A calibration control subroutine 140 includes programcode for controlling the calibration tool 108 to determine thecharacteristics of the laser beam generated by the system 100 and toprovide feedback data to calibrate the laser 102 and delivery systemoptics 104 for improved scanning accuracy and precision. The calibrationdata can also be used to adjust the ablation algorithm in the scanninglocation determination subroutine 134 to revise the calculation ofscanning spot locations and numbers.

[0205] The computer program 130 can be used to determine the pulselocations and numbers, and to calibrate the laser system 100 in realtime. In this way, the target profile for a patient can be entered justprior to performing the corneal sculpting procedure.

III. EXAMPLES

[0206] The following experimental examples are used to illustrate themethods of the present invention.

[0207] A. Least Square Fitting Approach Using Analytical Lens Equation

[0208] The examples were undertaken using a scanning laser system havinga Questek impulse excimer with a pre-objective galvo scanner. The pulserate was 70 Hz. The treatment time is dependent on the amount of tissueremoved per pulse (or tissue removed per second). For instance, a 6 mmspot will remove 9 times more material per pulse than a 2 mm spot. Tokeep the ablation rate constant, the 2 mm spot will need to increase therepetition rate by 9 as compared to the 6 mm spot. Decreasing the spotfurther will increase the number of pulses dramatically, as illustratedin the following table. Volume Spot Size [mm^(3]) 0.5 mm 4.7 × 10⁻⁹ 1.0mm 1.9 × 10⁻⁸ 2.0 mm 7.5 × 10⁻⁸ 6.0 mm 6.8.10⁻⁷

[0209] 1. Myopia

[0210]FIG. 9 is a graph comparing a fit using the method of the presentinvention to the theoretical lens equation. The number of basis, M, usedin the fit is 27. The fit is given a standard deviation, σ=0.25, for thefirst 150 points and a value four times larger for the last 150 points,σ=1.0. Having σ smaller for the inner portion of the curve gives thecentral 3 mm of a 6 mm ablation more importance than the outer portion.The difference in the quality of fit between the inner and the outerportions is evident, where the fit is much better in the inner portion.The outer portion shows a slight departure.

[0211] One gross feature observed in this fit is the “bump” at theone-millimeter distance. This occurs because of the size of the spot (2mm) is an integral multiple of the ablation zone (6 mm). One way to fixthe problem is to increase the ablation zone slightly such that theablation zone is not an integer multiple of the spot. This will alsotaper the outer perimeter giving a smooth transition from the ablatedarea to the non-ablated area. Further, there are numerous small bumpsnear the center ablation zone in FIG. 9. These manifest themselves asconcentric rings when the ablations are tested on plastics. Increasingthe number of basis functions will eliminate most of the rings.

[0212] Shown in FIG. 10 is the same fit but with an increase from 27 to47 basis functions. The ablation zone is not an integer multiple of thespot. The fit has the extended smooth transition between the ablated andnon-ablated zones. Most of the undulations have disappeared. It shouldbe noted, however, that care needs to be taken when increasing thenumber of basis functions, because too many basis function may cause thesolution to require negative values making such solutions not physicallyobtainable.

[0213]FIG. 11 is an interferometer image of a −4 diopter ablation on aplastic using the simulation output from FIG. 10 (extended fit with 47basis functions) and a spot size of approximately 2.1 mm. The salientfeature is the circular over-ablation at 1 mm radius.

[0214]FIG. 12 shows the characteristics of the simulated ablationprofiles of a −4 diopter ablation with 47 basis functions withvariations in spot sizes from a 2 mm reference spot size. The curvesabove the 2 mm reference curve have 5, 10, 15 percent larger spot sizes(i.e., 2.1 mm, 2.2 mm, and 2.3 mm), while the curves below have 5, 10,15 percent smaller spot sizes (i.e., 1.9 mm, 1.8 mm, and 1.7 mm). Thedip or peak at 1 mm becomes larger as the spot size deviates fromnominal, and the under- or over-ablation is non-linear. As a result, achange in spot size may remove less or more material and change theoverall profile.

[0215] 2. Hyperopiia

[0216]FIG. 13 is an interferometer image of a hyperopic ablation profileof a +4 diopter ablation with 33 basis functions using a 2 mm spot sizefor a 6 mm X 9.6 mm ablation zone. The method to generate hyperopiclenses is similar to that for myopia. The fitting requirement forhyperopia is less restrictive than for myopia. In myopia, both theablation's central and perimeter contours (and everything in between)desirably follow the lens equation. In hyperopia, the ablation is anegative lens, so that the central portion has very few pulses and theperimeter has the majority of the pulses. The fitting in the centralportion is much less stringent for hyperopia than for myopia. The pulsesdesirably are able to extend beyond the 6 mm treatment zone length inhyperopia. The outermost pulses have their centers at the 3 mm radius.For a 2 mm spot, it is desirable for the pulse to extend beyond the 3 mmradius going to about 3.5 mm to achieve a good fit. This extensionproduces an overall ablation diameter of about 9.6 mm. To reduce theablation size, a smaller spot may be used.

[0217]FIGS. 14A and 14B are graphs of two hyperopic lenses employing thesame basis functions used in the myopia examples and two different spotsizes. FIG. 14A uses a 2 mm spot and FIG. 14B uses a 1.5 mm spot. Theeffective dimensions are 6×10 and 6×8 mm ablations, respectively. Atransition zone is observed from the 3 mm radius onwards. The spot sizehas a large influence on this transition zone. The larger the spot size,the larger the transition zone becomes. For a 6 mm optical zone, a 2 mmspot produces a transition zone of 4.8 mm radius, while a 1.5 mm spotproduces a transition zone of 3.9 mm radius.

[0218] Fits to a hyperopic lens produce different characteristics fromthose obtained from the myopic fits. The hyperopic fits are relativelysimple because the perimeters are allowed to “float” thereby producing atransition zone. As with myopia, an extended zone is desirable for thehyperopic fit. This extended zone is about 0.8 times the radius. Thisplaces the spots at 3.8 mm for a 2 mm spot. Adding the radius of the 2mm spot gives the ablation zone radius of 4.8 mm. $\begin{matrix}{{treatment}\quad {zone}} \\{radius}\end{matrix} = {{{optical}\quad {zone}} + {1.8\quad \left( {{spot}\quad {{size}/2}} \right)}}$

[0219] Based on this equation, the treatment zone diameter for a 1.5 mmspot is 8.7 mm.

[0220] 3. Hyperopic Cylinder

[0221] Using the appropriate lens equation and basis functions, ahyperopic cylinder ablation profile is generated for a +4 diopterablation with 33 basis functions using a 2 mm spot size for a 6 mm×9.6mm ablation zone, as shown in FIG. 15.

[0222] Employing the method of the invention, the fit returns theamplitudes of the basis functions. Stacking the pulses vertically willgive the correct depth but it will not produce the lateral dimensionsneeded. As described earlier, offsetting the pulses will produce thetrench. The boundary conditions are given by the rectangular perimeterof the hyperopic cylinder. For example, a 5×9 mm hyperopic cylinder hasa 9 mm width along the lens direction (similar to plain hyperopia) andhas a 5 mm cut along the perpendicular directions. The 5 mm is thelength of the trench not including the transition zone. The transitionzone is the overlap-offset used to generate the basis “trenches.” Theslop of the transition zone is given by${slope} = \frac{{depth}\quad {per}\quad {pulse}}{offset}$

[0223] The offset is constrained by the total depth, the spot size, andtrench length. The relationship among the total depth, spot size, andoffset is given by${{step}\quad {size}} = \frac{{spot}\quad {size}}{{depth}\text{/}{depth}\quad {per}\quad {pulse}}$

[0224] The transition zone is one spot size in width. The number of thetotal pulses for this basis trench is not simply the depth/depth perpulse, since it also depends on the trench length. Instead, the numberof pulses is${{{no}.\quad {of}}\quad {pulses}} = {\frac{{length}\quad {of}\quad {trench}}{spotsize}*\frac{depth}{{depth}\quad {per}\quad {pulse}}}$

[0225] For myopia and hyperopia, the periodicity is reduced considerablyby adding a random starting phase to each ring. This approach is notavailable to ablation profiles with cylindrical symmetry such ashyperopic and myopic cylinders. Starting the trench with a repetitiousposition will produce a deeper pattern at that position with afeathering effect at the ending positions. The feathering effect occursbecause the number of pulses per basis is different making the offsetsdifferent. A very small random phase or position is desirably added toreduce the periodicity of the scans. For instance, a small random edgesposition is added to the locations to blend the spot.

[0226] Generating a myopic cylinder ablation profile is similar togenerating a hyperopic cylinder profile. The lens equations aredifferent but the same basis functions are used. Fitting the basisfunctions to the lens equation for a myopic cylinder is more stringentthan for a hyperopic cylinder, however, because the fit is preferablygood from the center to the perimeter with a small allowance for thetransition zone. The fit in the center region for a hyperopic cylinderis less rigorous.

[0227] In addition, care needs to be taken when the treatment ringapproaches the origin. Analytically, the basis functions for the myopiccylinder are equivalent to the ones used for hyperopic cylinder. Thedifference is the hyperopic cylinder treatment does not approach theorigin. At the origin, the basis functions overlap on themselves. Theoverlap is $\begin{matrix}{{y_{i}(x)} = {\sqrt{\left( {{spot}\quad {{size}/2}} \right)^{2} - \left( {x - x_{0i}} \right)^{2}} +}} \\{\sqrt{\left( {{spot}\quad {{size}/2}} \right)^{2} - \left( {x + x_{0i}} \right)^{2}}}\end{matrix}$

[0228] where spot size/2<(x_(0i)+x). This “overlap” is part of thealgorithm and is used if the fit goes from 0 to the treatment zoneradius, to account for the situation where the spot passes through theoriginal ([0]).

[0229] B. Simulated Annealing Approach

[0230] An example employs the simulated annealing approach for fittingthe ablation of a −2 diopter elliptical profile which has a major axislength of 16 mm and a minor axis length of 12 mm, and a maximum depth ofabout 31 μm. Due to the bilateral symmetry, only one quarter of thesurface is used in the fitting calculations. Any overlaps caused bylarge pulses near the axes of symmetry are reflected back relative tothe axes of symmetry to ensure that the calculations are correct. Thepulse shapes are circular having diameters of 5.5 mm, 5 mm, 4 mm, 3 mm,2 mm, and 1 mm at evenly spaced (100 μm) intervals. There are about 2000total pulse basis functions. The merit function in this case is the rootmean square error between the target shape and the current solution.

[0231]FIGS. 16a-16 c show the simulated annealing process at differentstages of the iteration to achieve the best fit. Each figure shows (1) aworking solution obtained from the current best guess of the treatmentplan or table, (2) the target profile, and (3) delta representing thedifference between the working solution and the target profile. At theinitial guess with zero iteration shown in FIG. 16a, the errors arequite large. The errors are reduced to about 2 μm or less after 10,000iterations as shown in FIG. 16b. At 1,000,000 iterations in FIG. 16c,the errors are less than 1 μm and well distributed over the surface. Thecontour plots of the target profile and the calculated profile areillustrated in FIGS. 17a and 17 b, respectively. About 50% of thecalculated profile falls within 0.1 μm of the target profile, and about80% of the calculated profile falls within 0.25 μm of the targetprofile. The solution is used by a Tencor scanner to produce theablation profile shown in FIG. 18.

[0232] It is to be understood that the above description is intended tobe illustrative and not restrictive. Many embodiments will be apparentto those of skill in the art upon reviewing the above description. Forinstance, some of the above equations are developed using a fixed spotshape and size. In alternate embodiments, the spot shape and/or spotsize may be varied. Variation of spot size and shape is described, forinstance, in commonly assigned European Patent No. 0628298 for “METHODAND SYSTEM FOR LASER TREATMENT OF REFRACTIVE ERRORS USING OFFSETIMAGING”, filed May 5, 1994 (Application No. 94303256.5) and issued Apr.1, 1998, which is incorporated herein by reference in its entirety.Different merit functions can be prescribed. Other numerical approachesmay be used to solve the target profile fitting equation. The desiredlens profile may be an elliptical profile, a hyperopic ellipticalprofile, a myopic elliptical profile, a circular profile, a linearprofile, an asymmetric profile, or an arbitrary two-dimensional lensprofile. The scope of the invention should, therefore, be determined notwith reference to the above description, but instead should bedetermined with reference to the appended claims along with their fullscope of equivalents.

What is claimed is:
 1. A method of generating a treatment table forablating tissue using a scanning laser beam for generating scanningspots over a treatment region larger in area than the scanning spots,the method comprising: providing a target function representing adesired lens profile for ablating the tissue by scanning spots of thelaser beam on the tissue; providing a basis function representing atreatment profile produced by scanning with overlapping scanning spotsof the laser beam din a treatment pattern; and fitting the targetfunction with the basis function to obtain a treatment table includingscanning spot locations and characteristics for the overlapping scanningspots of the laser beam.
 2. The method of claim 1 wherein the basisfunction is a two-dimensional function representing a two-dimensionalsection of a three-dimensional treatment profile which has symmetry withrespect to the two-dimensional section extending along the treatmentpattern.
 3. The method of claim 2 wherein the treatment pattern isgenerally linear or generally circular.
 4. The method of claim 1 whereinthe target function is a two-dimensional function representing atwo-dimensional section of a three-dimensional lens profile which hassymmetry with respect to the two-dimensional section extending along thetreatment pattern.
 5. The method of claim 4 wherein the target functionrepresents an ablation depth as a function of a distance from an opticalaxis of a cornea.
 6. The method of claim 1 wherein fitting the targetfunction with the basis function includes fitting at N discreteevaluation points.
 7. The method of claim 6 wherein the basis functionincludes M discrete basis functions representing M overlapping scanningspots.
 8. The method of claim 7 wherein the M discrete basis functionsrepresent M overlapping scanning spots across a treatment zone lengthrepresenting the length across a generally two-dimensional section whichis oriented normal across a generally straight treatment pattern orwhich is oriented radially across a generally circular treatmentpattern.
 9. The method of claim 8 wherein the scanning spots aregenerally circular and have a generally uniform energy profile.
 10. Themethod of claim 9 wherein (A) for a treatment profile having a generallyuniform two-dimensional section oriented normal across a generallystraight treatment pattern, the discrete basis functions represent thetwo-dimensional section as X _(i)(x _(j))=y _(i)(x _(j))={squareroot}{square root over ((s/2)²−(x _(j) −x _(0i))²)} or (B) for atreatment profile having a generally uniform two-dimensional sectionoriented radially across a generally circular treatment pattern, thediscrete basis functions represent the two-dimensional section as${X_{i}\left( x_{j} \right)} = {{\theta_{i}\left( x_{j} \right)} = {\cos^{- 1}\left( \frac{x_{j}^{2} + x_{0i}^{2} - \left( {s/2} \right)^{2}}{2 \cdot x_{0i} \cdot x_{j}} \right)}}$

 where s is the diameter of the scanning spot; j=1, . . . , N; x_(j) isa reference x-coordinate for the two-dimensional section measured froman optical axis of the cornea of a j^(th) evaluation point for thecenter of the scanning spot; x_(0i) is an x-coordinate for a center ofan i^(th) scanning spot; (x_(0i)−s/2)≦x_(j)≦(x_(0i)+s/2); y_(i)(x_(j))is a depth of the i^(th) basis function for the generally straighttreatment pattern; and θ_(i)(x_(j)) is a coverage angle of the i^(th)basis function for the generally circular treatment pattern.
 11. Themethod of claim 10 wherein x_(0i) is specified for M number of equallyspaced scanning spots as x_(0i)=i*[(L−s+e)/M], where L is the treatmentzone length; e is an extended zone; and i=1, . . . , M.
 12. The methodof claim 11 wherein e is set to about 0.1 to about 0.5 mm.
 13. Themethod of claim 7 wherein M is equal to about 7 to about
 97. 14. Themethod of claim 7 further comprising refitting the target function withthe basis function by varying the number of scanning spots M to iteratefor a best fit.
 15. The method of claim 6 wherein the target functionis: (A) for myopia and myopic cylinder, $\begin{matrix}{{f\left( x_{j} \right)} = {\sqrt{R_{1}^{2} - x_{j}^{2}} - \sqrt{\left( \frac{R_{1}\left( {n - 1} \right)}{n - 1 + {R_{1}D}} \right) - x_{j}^{2}} + C}} & {or}\end{matrix}$

(B) for hyperopia and hyperopic cylinder, $\begin{matrix}{{f\left( x_{j} \right)} = {R_{1} - \frac{R_{1}\left( {n - 1} \right)}{n - 1 + {R_{1}D}} - \sqrt{R_{1}^{2} - x_{j}^{2}} + \sqrt{\left( \frac{R_{1}\left( {n - 1} \right)}{n - 1 + {R_{1}D}} \right) - x_{j}^{2}}}} & {or}\end{matrix}$

(C) for phototherapeutic keratectomy, f(x _(j))=d;  where0≦x_(j)≦(L−shift); j=0, 1, . . . , N−1;${C = {\sqrt{R_{1}^{2} - {s^{2}/4}} + \sqrt{\left( \frac{R_{1}\left( {n - 1} \right)}{n - 1 + {R_{1}D}} \right) - \frac{s^{2}}{4}}}};$

x_(j) is an x-coordinate measured from an optical axis of the cornea ofthe j^(th) evaluation point for the center of the scanning spot; s isthe diameter of the scanning spot; R₁ is the anterior radius ofcurvature of the cornea in meters; R₂ is the final anterior radius ofcurvature of the cornea in meters; n=1.377 is the index of refraction ofthe cornea; D is the lens power of the scanning spot in diopters; L isthe treatment zone length representing the length across a generallyuniform section which is oriented normal across a generally straighttreatment pattern for myopic or hyperopic cylinders, or which isoriented radially across a generally circular treatment pattern formyopia or hyperopia; shift is the amount of emphasis shift; and d is aconstant depth.
 16. The method of claim 15 wherein the shift is about 0to about 0.2.
 17. The method of claim 15 wherein x_(j)=j*[(L−shift)/N].18. The method of claim 15 wherein the basis function includes Mdiscrete basis functions representing M overlapping scanning spots, andwherein fitting the target function with the basis function comprisessolving the following equation for coefficients a_(i) representingtreatment depth for the i^(th) scanning spot:${f\left( x_{j} \right)} = {\sum\limits_{i = 1}^{M}{a_{i}{X_{i}\left( x_{j} \right)}}}$

where X_(i)(x_(j)) is the i^(th) basis function; and i=1, . . . , M. 19.The method of claim 6 wherein fitting the target function and the basisfunction comprises specifying a deviation for each of the N discreteevaluation points.
 20. The method of claim 19 further comprisingrefitting the target function with the basis function by varying thedeviations to iterate for a best fit.
 21. The method of claim 1 whereinfitting the target function and the basis function comprises evaluatingcloseness of the fit and repeating the fitting step if the closenessdoes not fall within a target closeness.
 22. The method of claim 1wherein the target function and the basis function are fitted using aleast square fit.
 23. The method of claim 1 further comprisingrandomizing the scanning spot locations of the treatment table toproduce a random scanning order.
 24. The method of claim 1 furthercomprising refitting the target function with the basis function byvarying the size of at least one of the scanning spots to iterate for abest fit.
 25. The method of claim 1 wherein the scanning spotcharacteristics of a scanning spot at a scanning spot location includeshape, size, and depth of the scanning spot at the scanning location.26. The method of claim 1 wherein the scanning spots have differentsizes.
 27. The method of claim 1 further comprising specifying thetreatment pattern for scanning with overlapping scanning spots of thelaser beam.
 28. The method of claim 1 wherein the target function andthe basis function are fitted using a simulated annealing process. 29.The method of claim 1 further comprising specifying a merit functionrepresenting an error of fit between the target function and the basisfunction; and minimizing the merit function.
 30. The method of claim 1further comprising specifying a merit function representing an error offit between the target function and the basis function; monitoring atotal number of the scanning spots in the treatment table; andminimizing the merit function and the total number of the scanning spotsin the treatment table.
 31. The method of claim 1 further comprisingrefitting the target function with the basis function by selecting ascanning spot location and varying the characteristics of the scanningspot at the selected location to iterate for a best fit.
 32. A method ofgenerating a treatment table for ablating tissue using a scanning laserbeam for generating scanning spots over a treatment region larger inarea than the scanning spots, the method comprising: providing a lensfunction representing a desired lens profile for ablating the tissue byscanning spots of the laser beam on the tissue; providing a basisfunction representing a treatment profile produced by the overlappingscanning spots along a treatment path, the basis function representing asection oriented across the treatment path; and fitting the lensfunction with the basis function to obtain a treatment table includingscanning spot locations and characteristics for the overlapping scanningspots of the laser beam.
 33. The method of claim 32 wherein the scanningspots are generally circular and have a generally uniform energyprofile, and the basis function includes M discrete basis functionsrepresenting M overlapping scanning spots.
 34. The method of claim 33wherein the treatment profile is symmetrical with respect to an axis ofsymmetry, and the discrete basis functions are${\theta_{i}(x)} = {\cos^{- 1}\left( \frac{x^{2} + x_{0i}^{2} - \left( {s/2} \right)^{2}}{2 \cdot x_{0i} \cdot x} \right)}$

where s is the diameter of the scanning spot; x is an x-coordinatemeasured from the axis of symmetry; x_(0i) is an x-coordinate for acenter of an i^(th) scanning spot; (x_(0i)−s/2)≦x≦(x_(0i)+s/2); andθ_(i)(x) is a coverage angle of the i^(th) basis function.
 35. Themethod of claim 34 wherein x_(0i) is specified for M number of equallyspaced scanning spots as: x _(0i) =i*[(L−s+e)/M], where L is thetreatment zone length of the section oriented radially across thetreatment profile; e is an extended zone; and i=1, . . . , M.
 36. Themethod of claim 34 wherein fitting the lens function with the basisfunction comprises solving the following equation for coefficients a_(i)representing treatment depth for the i^(th) scanning spot:${f(x)} = {\sum\limits_{i = 1}^{M}{a_{i}{X_{i}(x)}}}$

where f(x) is the lens function; and i=1, . . . , M.
 37. The method ofclaim 36 wherein the lens function is: (A) for myopia,${f(x)} = {\sqrt{R_{1}^{2} - x^{2}} - \sqrt{\left( \frac{R_{1}\left( {n - 1} \right)}{n - 1 + {R_{1}D}} \right) - x^{2}} + {C\quad {or}}}$

(B) for hyperopia, $\begin{matrix}{{f(x)} = {R_{1} - \frac{R_{1}\left( {n - 1} \right)}{n - 1 + {R_{1}D}} - \sqrt{R_{1}^{2} - x^{2}} + \sqrt{\left( \frac{R_{1}\left( {n - 1} \right)}{n - 1 + {R_{1}D}} \right) - x^{2}}}} & {or}\end{matrix}$

(C) for phototherapeutic keratectomy, f(x)=d;  where 0≦x≦(L−shift);${C = {\sqrt{R_{1}^{2} - {s^{2}/4}} + \sqrt{\left( \frac{R_{1}\left( {n - 1} \right)}{n - 1 + {R_{1}D}} \right) - \frac{s^{2}}{4}}}};$

s is the diameter of the scanning spot; R₁ is the anterior radius ofcurvature of the cornea in meters; R₂ is the final anterior radius ofcurvature of the cornea in meters; n=1.377 is the index of refraction ofthe cornea; D is the lens power of the scanning spot in diopters; L isthe treatment zone length; shift is the amount of emphasis shift; and dis a constant depth.
 38. The method of claim 36 further comprisingdividing the depth (a_(i)) for the i^(th) scanning spot by a depth perpulse of the laser beam to obtain a number of pulses per an i^(th)treatment ring for the i^(th) scanning spot; and dividing the number ofpulses per treatment ring by 2π to obtain an angular spacing betweenpulses for the i^(th) treatment ring.
 39. The method of claim 32 whereinthe scanning spots have a fixed spot size and a fixed spot shape. 40.The method of claim 32 wherein at least one of the spot size and spotshape of the scanning spot is variable.
 41. A method for fitting athree-dimensional target profile, the method comprising: providing atwo-dimensional basis function including overlapping portions torepresent a three-dimensional profile which has symmetry with respect toa two-dimensional section extending along a treatment pattern; andfitting the three-dimensional target profile with the two-dimensionalbasis function to obtain a distribution of the overlapping portions. 42.The method of claim 41 wherein the three-dimensional profile hassymmetry with respect to a two-dimensional section oriented radiallyfrom an axis of symmetry and extending in a generally circular treatmentpattern around the axis.
 43. The method of claim 42 wherein theoverlapping portions are generally circular, and the two-dimensionalbasis function comprises discrete basis functions each representing acoverage angle of one of the overlapping portions as a function of adistance from the axis of symmetry.
 44. The method of claim 41 whereinthe three-dimensional profile has symmetry with respect to atwo-dimensional section oriented normal across a generally straighttreatment pattern.
 45. The method of claim 44 wherein the overlappingportions are generally circular, and the two-dimensional basis functioncomprises discrete basis functions each representing a depth of one ofthe overlapping portions as a function of a distance from the axis ofsymmetry.
 46. A system for ablating tissue, the system comprising: alaser for generating a laser beam; a delivery device for delivering thelaser beam to a tissue; a controller configured to control the laser andthe delivery device; and a memory, coupled to the controller, comprisinga computer-readable medium having a computer-readable program embodiedtherein for directing operation of the system, the computer-readableprogram including a first set of instructions for generating a treatmenttable including scanning spot locations and characteristics for ablatingthe tissue over a treatment region larger in area than the spot size ofthe laser beam to achieve a desired lens profile for ablating thetissue, a second set of instructions for controlling the laser togenerate the laser beam, and a third set of instructions for controllingthe delivery device to deliver the laser beam to the tissue according tothe treatment table.
 47. The system of claim 46 wherein the first set ofinstructions of the computer-readable program includes: a first subsetof instructions for providing a target function representing the desiredlens profile for ablating the tissue by scanning spots of the laser beamon the tissue; a second subset of instructions for providing a basisfunction representing a treatment profile produced by the overlappingscanning spots in a treatment pattern; and a third subset ofinstructions for fitting the target function with the basis function toobtain the treatment table including the scanning spot locations andcharacteristics for the overlapping scanning spots of the laser beam.48. The system of claim 47 wherein the second subset of instructionsprovide a basis function which is a two-dimensional functionrepresenting a two-dimensional section of a three-dimensional treatmentprofile having symmetry with respect to the two-dimensional sectionextending along the treatment pattern.
 49. The system of claim 47wherein the first set of instructions of the computer-readable programincludes a fourth subset of instructions for refitting the targetfunction with the basis function by varying the spot size of the laserbeam to iterate for a best fit.
 50. The system of claim 47 wherein thefirst set of instructions of the computer-readable program includes afifth subset of instructions for evaluating closeness of the fit andrepeating the fitting step if the closeness does not fall within atarget closeness.
 51. The system of claim 47 wherein the first set ofinstructions of the computer-readable program includes a sixth subset ofinstructions for randomizing the scanning spot locations for thetreatment table to produce a random scanning order.
 52. The system ofclaim 47 wherein the first set of instructions of the computer-readableprogram includes a seventh subset of specifying the treatment patternfor scanning with overlapping scanning spots of the laser beam;
 53. Thesystem of claim 47 wherein the scanning spot characteristics of ascanning spot at a scanning location include shape, size, and depth ofthe scanning spot at the scanning location.
 54. The system of claim 47wherein the desired lens profile is selected from the group consistingof an elliptical profile, a hyperopic elliptical profile, a myopicelliptical profile, a circular profile, and a linear profile.
 55. Thesystem of claim 47 wherein the desired lens profile is asymmetric. 56.The system of claim 47 wherein the desired lens profile comprises anarbitrary two-dimensional lens profile.